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ABSTRACT 

We measure the turbulent diffusion coefficient of dust grains embedded in 
magnetorotational turbulence in a protoplanetary disc directly from numerical 
simulations and compare it to the turbulent viscosity of the flow. The simulations 
are done in a local coordinate frame comoving with the gas in Keplerian rotation. 
Periodic boundary conditions are used in all directions, and vertical gravity is not 
applied to the gas. Using a two-fluid approach, small dust grains of various sizes 
(with friction times up to i7oTf = 0.02) are allowed to move under the influence 
of friction with the turbulent gas. We measure the turbulent diffusion coefficient 
of the dust grains by applying an external sinusoidal force field acting in the 
vertical direction on the dust component only. This concentrates the dust around 
the mid-plane of the disc, and an equilibrium distribution of the dust density is 
achieved when the vertical settling is counteracted by the turbulent diffusion away 
from the mid-plane. Comparing with analytical expressions for the equilibrium 
concentration we deduce the vertical turbulent diffusion coefficient. The vertical 
diffusion coefficient is found to be lower than the turbulent viscosity and to have 
an associated vertical diffusion Prandtl number of about 1.5. A similar radial 
force field also allows us to measure the radial turbulent diffusion coefficient. 
We find a radial diffusion Prandtl number of about 0.85 and also find that the 
radial turbulent diffusion coefficient is around 70% higher than the vertical. As 
most angular momentum transport happens through magnetic Maxwell stresses, 
both the vertical and the radial diffusion coefficients are found to be significantly 
higher than suggested by the angular momentum transport by Reynolds stresses 
alone. We also find evidence for trapping of dust grains of intermediate friction 
time in turbulent eddies. 
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1. INTRODUCTION 

Knowledge of the transport properties of particles embedded in a turbulent gas medium 
is important in many aspects of protoplanetary disc modeling. If the spatial number density 
distribution of dust grains in a disc is required for the model, one must know the effect of 
turbulent diffusion on the dust grains. 

Vertical diffusion — The distribution of tiny dust grains, with radii smaller than around 
100 fim, determines the observability of protoplanetary discs around young stellar objects 
through their contribution to the infrared parts of the spectrum. An interesting observational 
effect of turbulent diffusion is its influence on the vertical settling of dust grains. The 
settling affects the spectral energy distribution of protoplanetary discs, since flaring discs, 
i.e. where the scale height of the gas density increases with radial distance, have a much 
stronger mid- to far-infrared excess than self-shadowing discs, where the scale height after 
a certain distance from the protostar begins to fall with radial distance (e.g. DuUemond 
2002). Recent model calculations by Dullemond & Dominik (2004) show that the vertical 
settling of dust grains towards the mid-plane of the disc can change an initially flaring 
disc into a partially self-shadowing disc, thus effecting the observability of the disc. These 
calculations depend - among other things - on the strength of the turbulence in the disc (the 
turbulent viscosity) and on the turbulent diffusion coefficient of dust grains in the direction 
perpendicular to the disc mid-plane. Also, Ilgner et al. (2004) recently considered the effect 
of vertical mixing in protoplanetary discs on the distribution of various chemical species and 
found the distribution to be influenced greatly by mass transport processes, again underlining 
the importance of vertical turbulent diffusion in the modeling of protoplanetary discs. 

Radial diffusion — Crystalline silicate dust grain features observed in comet spectra 
are often attributed to radial mixing in the solar nebula (e.g. Hanner 1999). Silicate dust 
grains are formed primarily in amorphous form, but they can become crystalline if exposed 
to temperatures above ~ 1000 K. Such a heating can obviously have occurred in the hot 
inner parts of the solar nebula, whereas comets are expected to have formed in the cold outer 
regions of the nebula, so in this picture some radial mixing must take place between the inner 
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and outer nebula. Bockelee-Morvan et al. (2002) consider disc models where crystallization 
of silicates happens in the inner, hot parts of the disc. They calculate that in a few times 
10^ years the crystalline silicate fraction reaches a uniform value outside the crystallization 
region due to radial turbulent diffusion, and that the value can approach unity for realistic 
disc parameters. From high resolution observations of three protoplanetary discs, van Boekel 
et al. (2004) find that the inner 1-2 AU of these discs contain a higher crystalline silicate 
fraction than the outer 2-20 AU. This supports the theory that crystalline dust grains form 
in the hot inner disc and are subsequently transported to the outer disc by turbulent gas 
motion. 

The existence of chondrules (millimeter-sized solid inclusions found in primitive mete- 
orites, see e.g. Norton 2002) is believed to be the result of collisions and coagulation of small 
dust grains (Blum & Wurm 2000). The size distribution of chondrules may be explained by 
selective sorting in the turbulent solar nebula (Cuzzi et al. 2001). The first step in planet 
formation is the build-up of kilometer-sized rocky and icy planetesimals (in the planetesimal 
hypothesis of Safronov 1969), either from sticking or due to a gravitational instability in 
the vertically settled dust layer. In the latter case, the equilibrium scale height of the dust 
layer is determined by the turbulent diffusion coefficient of the dust grains in the vertical 
direction (Cuzzi et al. 1993). An alternative planet formation hypothesis, the gravitational 
instability hypothesis (see Boss 2003, and references therein), states that planets form as a 
direct gravitational instability in the gas of a protoplanetary disc. The ability of a disc to 
undergo gravitational instability depends on its density and temperature structure, which is 
again dependent on the distribution and thus the turbulent transport of tiny dust grains. 

It is a modern paradigm of protoplanetary discs that shear instabilities in the gas flow 
lead to turbulence, which is again responsible for such diverse effects as heating, angular 
momentum transport and diffusion. The actual turbulence is often parametrized in a single 
parameter, the turbulent viscosity (which can be non-dimensionalized into the a-value of 
Shakura & Sunyaev 1973). This single parameter determines both heating, angular mo- 
mentum transport and diffusion. Candidates for protoplanetary disc turbulence are many. 
Most pronounced linear instabilities are vertical convective instability (Lin & Papaloizou 
1980) and the magnetorotational instability (Balbus & Hawley 1991), although the former 
has proved to lead to inward rather than outward transport of angular momentum (Ryu & 
Goodman 1992). Other instabilities have been proposed, such as the baroclinic instability of 
Klahr & Bodenheimer (2003) which must be non-linear according to Klahr (2004), a linear 
Rossby wave instability (Li et al. 2000) and a linear stratorotational instability (Dubrulle 
et al. 2005; Shalybkov & Riidiger 2005). 

Today's most accepted source of turbulence is magnetorotational turbulence (MRI). 
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For completely ionized discs, the emergence of self-sustained turbulence through the linear 
magnetorotational shear instability seems inevitable, both in local shearing box simulations 
(Brandenburg et al. 1995; Hawley et al. 1995) and in global accretion disc simulations (Ar- 
mitage 1998; Arlt & Riidiger 2001). The application of the ideal MHD equations to proto- 
planetary discs is only justified where the ionization fraction is relatively high (e.g. Fromang 
et al. 2002; Semenov et al. 2004). This may be given in the hot and dust-free inner parts of 
the disc, as well as away from the mid-plane of the disc and at large radial distances where 
the ionization is determined by cosmic ray and high energy photon penetration. In proto- 
planetary discs this had lead to the concept of a magnetically dead zone near the mid-plane 
of the disc where the ionization fraction is too low to sustain MRI. Fleming & Stone (2003) 
consider local shearing box simulations with a vertically dependent ionization fraction and 
find that some turbulent stresses are induced in the dead zone by the surrounding MRI 
turbulence. Thus angular momentum can be transported even in regions of the disc that are 
not magnetorotationally unstable. 

It is often assumed that turbulent transport takes place as diffusion. For dust grains, 
the turbulent flux is assumed proportional to the gradient of the dust-to-gas ratio (Dubrulle 
et al. 1995). Such a prescription does not per se determine a certain value for the turbulent 
diffusion coefficient. Hence it is often parametrized to be a scalar that is equal to the 
turbulent viscosity of the gas for tiny grains but falls gradually for larger and larger grain 
sizes (Cuzzi et al. 1993; Schrapler & Henning 2004). One argument for setting the turbulent 
diffusion coefficient equal to the turbulent viscosity is that the radial velocity fluctuations 
are the base of both (non-magnetic) angular momentum transport and diffusion (Tennekes 
& Lumley 1972, p. 143). Such an approach is simple to use, but caution should be taken 
regarding its validity, both regarding the numerical value of the turbulent diffusion coefficient 
and regarding the isotropy that is implicitly assumed when making it a scalar. 

The validity of the whole diffusion description must also be addressed An obvious cause 
of concern is the presence of dust-trapping mechanisms in the turbulent gas flow. Gas 
turbulence and global pressure gradients, e.g. from vertical and radial stratification, are the 
cause of two important trapping mechanisms. Whenever the gas is pressure-supported and in 
force balance, the embedded dust grains feel an excess force in the opposite direction to the 
gas pressure gradient, since they can never be in the same force equilibrium without pressure 
support. The dust grains thus feel an acceleration relative to the gas. This has various 
effects, e.g. vertical settling of the dust layer towards the mid-plane of protoplanetary discs 
or inward radial drift of dust grains if there is an outwards decreasing gas pressure in the 
disc, a notorious problem in planet formation (Weidenschilling 1977). The dust grains reach 
a terminal velocity when the friction with the gas balances out the acceleration due to the 
missing pressure gradient. The terminal velocity of very small dust grains is proportional to 
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the friction time. In non-turbulent disc models only global pressure gradients are present, but 
in a turbulent disc local, fluctuating regions of high and low pressure are expected to occur. 
Then dust grains continuously move up the local pressure gradient, and this contributes 
to the random motion of the grains, which is responsible for turbulent diffusion. Magnetic 
pressure gradients actually give the same effect, as we will show analytically in Sect. 2.3. A 
local concentration of dust grains can not be described as diffusion, so one of the goals of 
this paper is to test the validity of the global diffusion picture in the presence of turbulent 
pressure gradient trapping. 

For larger dust grains, where the friction time becomes comparable to the orbital period 
of the disc, another important dust-trapping mechanism sets in. Stationary rotational struc- 
tures in the gas (e.g. anticyclones) are given by an equilibrium between the global Coriolis 
force from the rotating disc and the centrifugal force of the rotating structure. As they enter 
such a structure, large dust grains experience a slow acceleration, due to drag forces with 
the gas. Rotating initially with the gas, but much slower, the Coriolis force dominates over 
the centrifugal force, and the dust grains are sucked into the eddy. This vortex trapping was 
proposed by Barge & Sommeria (1995), and has since then been subject of much theoretical 
investigation (e.g. Chavanis 2000; Johansen et al. 2004). The conclusions are that vortices 
are extremely efficient at trapping dust grains, and this efficiency may even explain how 
gas planets are formed before the dispersion of the gas disc (Klahr & Bodenheimer 2005). 
Vortex trapping would seem to be potentially even more threatening to the global diffusion 
description than pressure gradient trapping. 

In this paper we measure the turbulent diffusion coefficient of dust grains directly from 
numerical simulations of three-dimensional magnetorotational turbulence. We treat the 
physics of protoplanetary discs in the shearing sheet approximation, in which a local co- 
ordinate frame corotating with the disc is considered. Dust is added as an extra fluid that 
interacts with the gas through a drag force. The turbulent diffusion coefficient is measured 
by exposing the dust fluid to an external force fleld and then comparing the resulting dust 
density with analytical expressions derived with a parametrized diffusion term. By compar- 
ing the measured value to the turbulent viscosity we examine whether the two are indeed 
equal as is often assumed. We speciflcally address the question of whether the diffusion 
coefficient is isotropic by measuring diffusion in both the vertical and the radial direction. 
Finally we quantify the effect of dust-trapping mechanisms on the whole diffusion picture 
by examining correlations between turbulent gas features and the dust-to-gas ratio. 

The paper is built up as follows. In Sect. 2 we describe the dynamical equations for the 
motion of gas and dust and the computer code that we use to solve them numerically. Then 
we go into details in Sect. 3 about how we deduce the turbulent diffusion coefficient from 
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computer simulations by comparing the equilibrium dust density with analytical expressions. 
In Sect. 4 we describe the units and the boundary conditions of the simulations. The results 
are described in the following two sections; Sect. 5 describes the turbulent evolution of the gas 
while Sect. 6 describes the evolution of the dust, especially the measured turbulent diffusion 
coefficients and diffusion Prandtl numbers and the presence of dust-trapping mechanisms in 
the gas. Finally conclusions, discussions and some outlook to potential further investigations 
into the subject of turbulent diffusion of dust grains appear in Sect. 7. 



2. DYNAMICAL EQUATIONS 

In this section we present the dynamical equations we use for gas velocity, gas density, 
magnetic vector potential, dust velocity and dust density. We integrate the dynamical equa- 
tions using the Pencil Code^. This is a finite difference code that uses sixth order centered 
spatial derivatives and a third order Runge-Kutta time-stepping scheme. See Brandenburg 
(2003) for details on the numerical schemes and test runs. The Pencil Code solves the dynam- 
ical equations in their non- conservative form and gives very similar results to the ZEUS code 
for the statistical properties of MRI turbulence (see Balbus & Hawley 1998, and references 
therein) . 

The Pencil Code requires artificial diffusivity terms in the dynamical equations to stabi- 
lize the finite difference numerical scheme. For the purpose of calculating turbulent diffusion 
coefficients it is vital that we can reduce the artificial mass diffusion as much as possible 
in order to distinguish the measured turbulent diffusion from the imposed artificial diffu- 
sion. The biggest contribution to the turbulent transport of dust grains is expected to 
come from the fast and far moving large scales, so keeping the large scales unaffected by 
diffusion is important. To achieve this we use hyperdiffusivity terms in all the dynamical 
equations. Hyperdiffusivity involves replacing the regular diffusivity terms (involving sec- 
ond order derivatives) with differential operators that use higher order derivatives. This 
quenches unstable modes at the smallest scales of the simulation, while at the same time the 
large scales are kept unaffected by diffusivity. We have checked, by varying the value of the 
artificial diffusion coefficient, that hyperdiffusion does not have any effect on the turbulent 
diffusion coefficients that we measure. The use of hyperdiffusivity is discussed further in Ap- 
pendix A. There the hyperversions that we adopt for viscosity, mass diffusion and resistivity 
are also presented. 



^The code is available at 
http:/ /www. nordita.dk/software/pencil-code 
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In this section we also develop a method for being able to treat numerically the motion 
of very tiny dust grains with friction times much shorter than the computational time-step 
of the gas. This so-called short friction time approximation is presented and discussed in 
the last part of the section. 



2.1. Gas Dynamics 

We consider the motion of gas and dust in the shearing sheet approximation (e.g. Gol- 
dreich & Tremaine 1978; Brandenburg et al. 1995). Here a local coordinate frame corotating 
with the disc at a distance Tq from the central source of gravity is considered. The coordinate 
axes are defined as following. The x-axis points always away from the central gravity source, 
and the y-axis points in the direction of the Keplerian flow (as seen from a non-comoving 
frame). The 2;-axis points perpendicular to the disc along the direction of the angular ve- 
locity vector i7o of the orbital motion. In this frame, the Keplerian flow velocity field has 
the linearized form Uq = —^QQxy = u^y^'y. We choose to measure velocities relative to 
the Keplerian flow, u — uq ^ u. Such a transformation introduces new shear terms in the 
equation of motion, but it has the advantage that the Keplerian velocity is zero everywhere. 
The equation of motion relative to the main shear flow in the shearing sheet approximation 
is 

(jfl (ill 1 1 

— = -(u ■ V)u - 4°)— + f{u) --VP + -JXB + f^{u, p) . (1) 
ot ^ oy p p 

The first term on the right hand side of equation (1) is the advection due to any velocity 
relative to the main shear flow, while the second term covers the advection due to the shear 
flow. The function f{u) is defined as 

/ 2f]oUy \ 

f{u) = I -i|2oM, j (2) 

and is an effect of Coriolis force. The last three terms in equation (1) are the pressure gradient 
force, the magnetic Lorentz force (where the volume current density J is defined through 
Ampere's law V x S = poJ), and a hyperviscosity term based on the function defined in 
equation (A4). We ignore the effect of vertical gravity on the gas, because we are interested 
in the ideal case to measure the isotropy/non-isotropy of magnetorotational turbulence and 
the local transport properties of the gas without introducing additional isotropy-breaking 
effects. Ignoring the stratification effectively means that we are considering the disc close to 
the mid-plane where the vertical gravity is vanishing. Future work on the properties of dust 
diffusion in local shearing box simulations should take the vertical stratification of the disc 
into account. 
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The evolution of the gas density p is determined by the continuity equation 

| = -r|-pV-u-n.Vp + /o(p), (3) 

where the first term on the right hand side is again an effect of advection due to the main 
shear flow. The two next terms come from the standard term V ■ (pit) from the continuity 
equation. In the last term we include artificial mass diffusion through the function /d defined 
in equation (A6). An isothermal equation of state is used where the pressure depends on 
the density as P = clp. Here Cg is the constant sound speed. 

The induction equation determines the evolution of the magnetic vector potential A. 
Evolving the vector potential has the advantage over evolving the magnetic field S = V x A 
in that it maintains a solenoidal magnetic field (i.e. V • S = 0) at all times. The induction 
equation in the shearing sheet approximation is 

^ = + + + f.^iA) . (4) 

The first term on the right hand side of equation (4) is the advection due to the main 
shear flow, while the second is the so-called magnetic stretching term, another effect of shear 
(Brandenburg et al. 1995). The two last terms are the standard electromotive force and a 
resistivity term based on the function defined in equation (A8). 



2.2. Dust Dynamics 

The dust grains are considered as a fluid without any pressure support. Any pressure 
gradient force on the dust due to collisions between dust grains and gas molecules can also 
be ignored since the solid density of dust grains is so large that the resulting acceleration is 
negligibly small. 

In the fluid approach, the equation of motion for the dust velocity relative to the Kep- 
lerian flow is 

_ = -[w ■ V)w - u^y^^ + f{w) + f^{w,n) - —{w - u) +g{x,y,z) . (5) 

The first four terms on the right hand side appear similar here as in the gas momentum 
equation. The last two terms in equation (5) are the drag force and an externally imposed 
force field g that we use to drive a non-zero diffusion equilibrium in the dust density. This 
is explained in more detail in Sect. 3. 
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We let the dust and the gas couple through a drag force proportional to, but in the 
opposite direction of, the velocity difference between the dust and the gas. The strength of 
the drag force is characterized by the friction time Tf . Any relative motion between dust and 
gas is damped by the drag force with an e-folding time of Tf . The physics of the dust grain 
and the gas enters in the expression of the friction time. When the mean free path of the 
gas molecules is larger than the dust grain radius, the dust grain is in the Epstein regime 
(e.g. Weidenschilling 1977). Here the friction time of a spherical dust grain with radius a, 
and solid density p, can be expressed as 

= , (6) 

CsP 

where Cg is the sound speed in the surrounding gas and p is the gas density. The sound speed 
and the gas density are approximately constant in the unstratified and isothermal case, so 
we can treat Tf as constant that depends only on the given particle radius and solid density. 

Treating dust as a fluid is justified as long as the mean free path of the fluid constituents 
is smaller than the typical dimensions of the system. In the case of the gas, one compares the 
mean free path of the molecules with the thickness of the disc Hq. For the dust grains the 
collisions among the grains themselves are unimportant for determining a mean free path. 
Here the collisions with the gas molecules are the important effect. The mean free path for 
the dust grains can be defined as the distance one grain has to float with respect to the 
gas before it has lost a signiflcant fraction of its momentum. For a spherical grain moving 
with a speed w relative to the gas, this value can be determined as £ = WTf. The condition 
for treating dust as a fluid is then that i <^ Hq. Because all motions are subsonic, we can 
replace w by Cg as an upper limit and get the expression i^o^f "C 1 for the validity of the 
fluid approach. 

In a typical solar nebula type protoplanetary disc, the scale height is of the order of 
Ho ~ 10^^ cm at ro = 5 AU, while the gas density can be taken to po ~ 10~^° g cm~^ at the 
same radial distance. Then the connection between grain radius and dimensionless friction 
time is 

a, = f2QT{Ho— ~ lO^f^oTfcm. (7) 
p. 

This means that, as a rule of thumb, the value of the dimensionless friction time corresponds 
to the radius of the dust grain measured in meters. 

To preserve momentum the gas should be affected by a drag force /drag = ^"^i^ Pa/ p{u — 
w) from the dust. Here pd/p is the dust-to-gas ratio. We shall ignore the back-reaction drag 
force from the dust on the gas, because the dust-to-gas ratio is small in the early stages of 
a protoplanetary disc. 
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We represent dust mass density by the number density n of dust grains. The continuity 
equation for the dust number density n is 



where we use artificial diffusion, through the function fnin) defined in equation (A6), only 
to stabilize the numerical scheme. By varying the value of the artificial diffusion coefficient 
-D3, which is defined in Appendix A, we have made sure that adding artificial diffusion has no 
effect on the measured turbulent diffusion coefficients. The value of needed to stabilize 
the runs are for all runs several orders of magnitude below the measured turbulent diffusion 
coefficient. 

We now have two possibilities to solve the dust equation of motion (eq. [5]): For large 
particles, with friction times larger than the Courant time-step (^^o^f > 0.001, see Table 1), 
we can use the explicit integration scheme from the Pencil Code. But for the smaller particle 
cases {f2oT{ <^ 0.001), where the friction time is much shorter than the Courant time-step, 
we will use the short friction time approximation, a semianalytical time integration that is 
presented below. 



The radii of dust grains observed in protoplanetary discs are often on the order of 
micrometers or even nanometers. The friction time of microscopic dust grains in a proto- 
planetary disc is very short compared to the orbital period, around a few minutes for the 
location of Jupiter in a typical solar nebula. That is of course not a problem for nature, 
but the smallest scales of computer simulations are many orders of magnitude larger than in 
nature, and thus the computational time-step for an explicit code such as the Pencil Code 
is much larger than the friction time of tiny dust grains. This causes a potential problem in 
resolving both timescales at the same time. To follow the motion of the tiniest dust grains, 
applying the explicit integration scheme as used in the Pencil Code, a time-step must be 
chosen that is at least an order of magnitude shorter than the friction time. Thus the compu- 
tation time for simultaneously following the evolution of gas and dust becomes prohibitively 
long. One can now either use an implicit integration scheme, which would introduce further 
problems, and also make major changes in the code necessary, or one can use a kind of 
semianalytical integration scheme that works as follows. 

For very short friction times, the dust is able to settle to an equilibrium velocity, where 
the drag force is exactly balanced by the other force terms that are present in the dust 





2.3. Short Friction Time Approximation 
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equation of motion (eq. [5]), on a timescale that is much shorter than the computational 
time-step of the gas. Thus it is possible, under a few reasonable assumptions, to solve 
algebraically for the terminal dust velocity as a function of gas velocity and density. To do 
this, we first subtract the gas equation of motion (eq. [1]) from the dust equation of motion 
(eq. [5]). This results in an equation for the evolution of relative velocity w — u, 

— ^+ It;. V iu- u-V M+uf ^ = f{w-u) {w-u)+g+-{VP-JxB) , 

at ^ ay ti p 

(9) 

where the viscosity terms have been ignored since any real physical viscosity is expected 
to be orders of magnitude weaker than the other force terms. We now assume that the 
computational time-step of the gas is much longer than the friction time, 6t ^ Tf. Here 6t 
will be given by the Courant criterion. This criterion determines the maximum time-step 
that can be taken by an explicit numerical scheme without becoming unstable. The allowed 
time-step gets shorter with increasing grid resolution. With the condition 6t ^ Tf, all terms 
from the gas equation of motion can be considered to be constant for the duration of the 
acceleration of the dust grain to its terminal velocity. This specifically also applies to the 
pressure gradient force and the Lorentz force that are present also in equation (9). Then 
we can search for a time-independent equilibrium solution for w — u. We expect any time- 
independent solution of equation (9) to have a dust velocity that is very close to the gas 
velocity, because the short friction time couples the dust velocity strongly to the gas velocity. 
Setting therefore w = uin all other terms than the drag force term (this is legitimized below) 
leaves the algebraic equilibrium equation 

= --(yj -u) +g + -{VP - J X B) . (10) 

Tf p 



Solving for w then yields 

W = U + T{ 



g + -(VP - J X B) 

P 
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Reinserting this solution into equation (9) shows that it was reasonable to ignore all ad- 
vection, shear and Coriolis terms, while keeping the friction, gravity, pressure gradient and 
Lorentz terms, as long as the friction time is sufficiently short. 

This is the short friction time approximation. The specific form of the short friction 
time dust velocity approximation depends on the forces that are assumed to work on the 
gas and on the dust, so equation (11) is only valid for the specific choice of force terms 
that are considered in this work. The presence of gravity in the dust velocity approximation 
comes from only considering gravity to work on the dust. This is is good for the purpose of 
measuring the turbulent diffusion coefficient, whereas in nature gravity of course affects both 
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dust and gas. The gravity term would then drop out of the short friction time approximation, 
but it would reappear in the form of the vertical pressure gradient of the stratified gas. One 
must also take into consideration that the dust velocity in equation (11) is expressed as a 
function of the resolved part of the gas velocity only. All unresolved small scales would also 
contribute to the random motion of the tiny dust grains (as would Brownian motion), but 
the important scales for turbulent transport are the largest scales in the box, since they 
contribute most to the total gas velocity field. 

The equilibrium dust velocity given by the short friction time approximation ensures 
that the relative velocity between dust and gas does not change on timescales shorter than 
the computational time-step. That means that if the gas is being accelerated, then the same 
amount of acceleration must be working on the dust, and so the relative velocity between 
the dust and the gas stays constant until sufficient time has passed for the pressure gradient 
force and the Lorentz force to change at the computational timescale. 

The ith component of the Lorentz force appearing in equation (11) can be rewritten as 



where the first term in the parenthesis on the right hand side is due to magnetic pressure 
and the second to magnetic tension. This allows the ith component of the short friction time 
approximation dust velocity to be rewritten as 



Thus dust grains move relative to the gas not only because of (additional) gravity and 
(missing) pressure gradient force, but also due to (missing) magnetic pressure gradient force 
and (missing) magnetic tension. We shall still refer to the effect as pressure gradient trapping, 
even though the magnetic tension term in equation (13) does not mimic a pressure gradient. 

We must stress again that the short friction time approximation is only valid for small 
particles. If one considers larger bodies (e.g. > 1 mm at tq = 5 AU in a typical solar nebula), 
first the Coriolis forces and then the advective transport terms can no longer be ignored. For 
these particles we directly integrate the dust equation of motion (eq. [5]) together with the 
other dynamical equations. With even larger objects finally the fiuid approach fails as soon 
as i > Hq. In this case one has to apply a particle algorithm to follow the dust evolution 
(e.g. Klahr & Henning 1997). 




(12) 




(13) 
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3. 



DIFFUSION COEFFICIENT 



In this section we describe how we calculate the diffusion coefficient of dust grains 
embedded in a turbulent gas. We do this by comparing the results of numerical simulations 
with analytical solutions to the non-turbulent flow equations that include a parametrized 
diffusion. 

If the turbulent motion of the gas and the dust has not been resolved, the continuity 
equation of the dust would have to incorporate an explicit diffusion term, 



The gas flow is here assumed to be completely stationary, and the only effect of the non- 
resolved turbulence is through the parametrized diffusion term. The continuity equation is 
written in a conservative form where the diffusion flux is proportional to and in the opposite 
direction of the gradient of the dust-to-gas ratio. This is the way turbulent diffusion is 
normally assumed to act (see e.g. DubruUe et al. 1995). 

The task now is to find a way to extract Dt from the non-stationary turbulent motion 
found in computer simulations. This is only possible if Vra is not zero everywhere, as 
otherwise the diffusion coefficient does not enter equation (14) at all for a constant p. One 
can now either follow the time dependent diffusion of an initial dust concentration somewhere 
in the center of the box or look for a time independent equilibrium solution. The first 
approach has the disadvantage that it is difficult to obtain good statistics, as one has always 
a very special distribution with a distinct wavelength, whereas the turbulence could act on 
all scales. Therefore we use the latter possibility and search for an equilibrium solution where 
we can achieve much better statistics. 

We force an equilibrium solution with a non-zero dust density gradient by exposing 
the grains to an external force field g. Depending on its specific form, this force field 
will eventually result in an equilibrium where the pile-up of dust grains imposed by g is 
balanced completely by mass diffusion in the opposite direction. By comparing the analytical 
expression for the equilibrium dust number density, whose only free parameter is -Dt, to the 
equilibrium density obtained when the turbulence is resolved in computer simulations, we 
can derive the turbulent diffusion coefficient. We will often refer to the external force field 
simply as gravity because of the qualitative similarities to real gravity. 

First the diffusion in the z-direction is considered. Here we define a vertical gravity field 




(14) 




(15) 
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where kz = 2tt/Lz in order to have periodic boundaries in the vertical direction. Here Lz 
is the vertical extent of the box, and z is defined to lie in the interval between —\Lz and 
\Lz. Using periodic boundary conditions demands that we use a periodic force field in order 
to have a periodic equilibrium solution. The gravity field defined in equation (15) is linear 
around the mid-plane, as the gravity field normally considered for thin discs also is, but away 
from the mid-plane it becomes zero again on the top and bottom boundaries of the box. Such 
a force gives a periodic dust distribution to determine the turbulent viscosity coefficient from 
(we will show below that the equilibrium logarithmic dust density becomes cosinusoidal with 
z). For a normal thin disc vertical gravity field, Qz = —il^z, the equilibrium logarithmic 
dust density becomes quadratic with z, which is obviously not periodic. 

To find the equilibrium dust number density, we solve now equations (5) and (14) for 
dw/dt = dn/dt = u = = Wy = 0, Wz = Wz{z) and n = n{z). This yields the differential 
equation system 

= -Wz^^ - —Wz- gosm{kzz) , (16) 

OZ Tf 

where we neglect the p-dependence in the diffusion term, because the turbulent gas density 
fluctuations are very small, as expected in subsonic turbulence. For any sufficiently short 
friction time, the advection term in equation (16) can be safely ignored, leaving only the 
algebraic equation 

= Wz - go sm{kzz) (18) 

Tf 

with the solution 

Wz = -Tigosm{kzz) . (19) 

Inserting equation (19) into equation (16) shows that the advection term is fully negligible 
for Tfgokz < 1. 

The equilibrium solution to the continuity equation must now be able to continuously 
replace material that is being transported towards the mid-plane by new material transported 
away from the mid-plane by diffusion. It is seen that equation (17) has the general solution 

lnn = ^ /" Wz{z)dz (20) 
Dz J 

for any integrable function Wz{z). Here we have assumed that there is no net flux of dust 
grains = by setting the contents of the parenthesis on the right hand side of equation 
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(17) equal to zero. Inserting the equilibrium dust velocity from equation (19) into the integral 
in equation (20) gives the equilibrium logarithmic dust number density as 

lnn = lnniH — cos(kzz) , (21) 

where In ni is an integration constant that corresponds physically to the logarithmic number 
density aX z = rL^L^. The amplitude of the cosine distribution depends only on friction 
time, gravity strength and gravity wave number, which are all known input parameters, and 
the unknown value of the turbulent diffusion coefficient in the vertical direction. Thus the 
value of the turbulent diffusion coefficient can be determined uniquely from this amplitude. 

The number density distribution in equation (21) is not normalized. The connection 
between rii and the column density Eq is 

^1/2 



rii exp 

1/2 



(t) 



27mi 1 
az = — / exp 



Tl9o /, X 

COS[KzZ) 



(t) 



d{kzz), (22) 



where the last equality holds because the cosine function is symmetric in z. The modified 
Bessel function of the first kind of order m is defined as 

J^(x) = i/ e"™^^cos(m0)d^, (23) 
Jo 

so the connection between and rii becomes simply 



Isolating rii finally yields 



ni = r- . (25) 



27r/o 



For infinite diffusion D^z'^ oo, the argument of the Bessel function in equation (25) is zero, 
and using /o(0) = 1 from equation (23), we get Eq = riiLz- Thus rii = uq, where no is the 
average dust number density in the box, as expected for the special case of infinite diffusion. 
In the case of a finite diffusion coefficient, ni 7^ no. 

For the radial x-direction, a similar sinusoidal gravity field can be defined to give the 
equilibrium dust density as 

In n = In ni H — cos^k^x) , (26) 
kxDx 

formally identical to the vertical case. The derivations are given in Appendix B. With equa- 
tions (21) and (26) we are armed with two powerful analytical expressions for the number 
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density distribution of dust grains in diffusion equilibrium with an externally imposed force 
field. By comparing computer simulations of magnetorotational turbulence with these ana- 
lytical results, we can extract the turbulent diffusion coefficient of the dust grains in both 
the vertical and the radial directions independently. The next sections describe the setup of 
the simulations and the results that we get. 

4. UNITS AND BOUNDARY CONDITIONS 

We adopt non-dimensional variables by measuring velocities relative to the isothermal 
sound speed, [u] = [w] = c^, and densities relative to the initial density in the box, [p] = p^; 
[n] = uq. The unit of dust-to-gas ratio ed is [ed] = eo = niQno/po, where nio is the mass of 
the individual dust grains. Time is measured in units of inverse Keplerian angular speed, 
[t] = f^Q^, although often stated in orbits T = 27ri7(^^. The unit of magnetic field is [B] = 
Cs y/poPo- Derived from these basic units are the unit of distance [x] = Cgf^Q^ and the unit of 
magnetic vector potential [A] = c^Q^^ ^ p^pQ. The unit of turbulent viscosity and turbulent 
diffusion coefficient can also be derived from the basic units to be \ut\ = [Dt] = c^f^Q^. In 
these units the turbulent viscosity and the turbulent a- value take the same numerical value. 

We choose a box length of 27r in all directions. In order to keep the background shear 
flow subsonic at all points we choose the arbitrary normalization = 0.2. We have checked 
by setting to unity that the evolution of the simulations indeed scale with the value of 
Hq, and thus that the scale- free diffusion coefficients and a- values are independent of the 
choice of i7o- 

Periodic boundary conditions are applied in all directions. Connected points at the 
periodic x-boundary have a time-dependent shift as is appropriate in the shearing sheet 
approximation. 

5. EVOLUTION OF GAS 

As an initial condition, we perturb the gas velocity components with random fluctu- 
ations of amplitude du ~ 10"^. The toroidal component of the magnetic vector potential 
is perturbed by a standing cosine wave Ay = Aq cos{k^x) cos{kyy) cos^k^z) of amplitude 
Aq = 0.2 and wave numbers kx = ky = kz = 1. The resulting vertical component of the 
magnetic field is = —AQkxSm{kxx) cos{kyy) cos^k^z) = Bo{x,y) cos^k^z). Such a wave 
is unstable to shear if k^ is sufficiently small (i.e. at sufficiently large scales). As shown by 
Balbus & Hawley (1991), the wave number interval for instability of the vertical magnetic 
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field component is < A;^ < \/?)Qq/vx, where the Alfven speed is defined as v\ = Bq/{^qPq). 
For < \Bq\ < 0.2, the upper hmit wave number is always larger than around = V^, so 
kz = 1, the largest scale present in the simulation, is well within the unstable regime. 

We run simulations in two different resolutions, 64^ and 128^. In Table 1 the parameters 
that are used in the different runs are listed. As there is no back-coupling from the dust 
on the gas, the gas evolution depends mainly on resolution, since the high resolution runs 
require less artificial diffusivity. The dust only affects the gas through its contribution to the 
computational time-step. 

5.1. Self-Sustained Turbulence 

Initially the magnetic and kinetic energies in the box increase, but the increase stops 
after around half an orbit, and then the magnetic and kinetic energies fall slowly to an 
equilibrium state during a transition time of around ten orbits. In the equilibrium state the 
turbulence is self-sustained, in the sense that energy is pumped from the gravitational field 
primarily into magnetic energy (via the magnetorotational instability). The Lorentz force 
transfers some of this magnetic energy into turbulent kinetic energy which again transfers 
energy back into the magnetic field in a dynamo process. Finally the energy is dissipated 
through resistivity and viscosity. The whole process is sketched in Brandenburg et al. (1995). 
Because we assume an isothermal equation of state, there is no heating of the gas due to 
dissipative processes. 

The time evolution of kinetic energy components, magnetic energy components, and 
Reynolds and Maxwell stresses is shown in Fig. 1 for a time span of 100 orbits. All turbulence 
parameters are approximately constant in time, within a certain fiuctuation interval, and 
show no sign of decaying after the steady state has set in after around ten orbits. Most of the 
kinetic energy (top panels) is present in the horizontal components of the velocity field, which 
is always measured relative to the Keplerian fiow, whereas the vertical component contains 
a factor of two lower kinetic energy (this anisotropic trend is normal to MRI simulations, 
see e.g. Hawley et al. 1995). For the magnetic energy (middle panels), almost the entire 
energy is kept in the toroidal component of the magnetic field. The ratio between kinetic 
and magnetic energies stays approximately constant in time with the magnetic energy being 
a factor of around two higher than the kinetic energy. The Reynolds and Maxwell stresses 
(shown in the two bottom panels) can be converted into a turbulent viscosity and normalized 
to a turbulent a- value of Shakura & Sunyaev (1973). These values are shown for the different 
runs in the second and third columns of Table 2. The magnetic a- value is around a factor 
of four times the non-magnetic, so most angular momentum transport happens because of 



- 18 - 



magnetic fields. In the shearing sheet approximation the Keplerian background velocity is 
linear in space, so there is no pile-up of angular momentum anywhere in the box. 

The magnetorotational instability injects energy at the largest scales of the box. The 
smaller scales are then set in motion as the large scale motion cascades down to smaller 
and smaller scales. Under the assumption that there is no pile-up of kinetic energy at any 
scales, the Fourier spectrum should obey a Kolmogorov-law u{k) ~ k~^^^. The Fourier 
spectra of all velocity components for 64"^ and 128'^ runs are shown in Figs. 2 and 3. For 
reference a k~^^^ line is shown. The spectra are averages taken from 10 to 100 orbits. At 
large scales, the power spectra approximately obey a Kolmogorov law, but at smaller scales, 
where dissipation becomes important, the slope becomes steeper. There is some excess power 
at the very smallest scales, especially for the 64^ run. This is due to unstable modes at the 
smallest scales of the box. Curiously the excess power is only present in the radial and 
vertical directions and not in the toroidal direction, but this may be an effect of the shearing 
out of all variables along y. The power in the small scale modes is still negligible compared 
to the large scales, so the rise in power does not influence the diffusion of the dust. These 
rises in power are typical for simulations with a low diffusivity, see e.g. Haugen et al. (2004). 
According to mixing length theory, the contribution from the different length scales to the 
total turbulent diffusion coefficient scales as ~ Uk/k ~ k~^^^, so the largest scales of 
the box are expected to give the dominating contribution to the total turbulent diffusion 
coefficient. 

One also sees from Figs. 2 and 3 that in both cases the vertical velocity amplitude on the 
large scales is smaller than the radial and toroidal velocity amplitude at large scales. This 
gives already a hint that vertical turbulent diffusion might be weaker than radial turbulent 
diffusion. 



6. EVOLUTION OF DUST 

The dust is initially at rest and has a constant number density n{x, y, z) = Uq. It is 
then set free to evolve under the influence of friction with the gas and the imposed gravity 
field. The dust begins to concentrate near the center of gravity (horizontal mid-plane, 
with z = 0, for vertical gravity, vertical mid-plane, with x = 0, for radial gravity), but 
eventually an equilibrium conflguration is reached where the turbulent diffusion prevents 
further concentration. This situation is shown in Fig. 4 for a 128"^ run with OqT{ = 2 x 10"'^ 
and vertical gravity. The run is labeled 128a_z in Table 1, and the friction time corresponds 
to tiny dust grains or molecules with radii of 0.2 micrometers in a typical solar nebula. The 
plot shows dust density contours at the sides of the simulation box. The turbulent motion 
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is clearly visible, and the resulting turbulent diffusion is the only reason why there is no 
further settling of the dust layer towards the mid-plane. The amplitude of the concentration 
around the mid-plane is maintained approximately constant for the entire duration of the 
simulation (one hundred orbits). 



6.1. Diffusion Timescale 

Before proceeding with measuring diffusion coefficients, we will first make an estimate of 
the time it takes to get from a constant dust density to the equilibrium where sedimentation 
is balanced by turbulent diffusion. We consider the case of vertical gravity. The logarithmic 
dust density must rise from an initial value In uq to the equilibrium value given by equation 
(21). When the amplitude of the equilibrium cosine function is small, Ainn <^ 1, then we 
can assume that ni ^ riQ. The increase in logarithmic dust density is then simply 

A Inn = cos{kzz) . (27) 

This increase is caused by the vertical sedimentation. In the short friction time approxima- 
tion the dust velocity can be written as Wz = —tiQq sm{kzz). The change in logarithmic dust 
density due to vertical settling can be approximated with the expression 

dlnn dwz 

= — ^ = kzTfgo cos{kzZ) . (28) 

Here we have ignored the advection of mass for simplicity. The diffusion timescale can 
now be estimated by dividing equation (27) with equation (28). This yields 

Rewriting the diffusion coefficient in dimensionless units as D^z^ = 6^z^c^f2Q^, the diffusion 
timescale can be written as f^o^D = [{cs^o^kz)'^S^^]~^ ■ With i7o = 0.2, kz = 1 and 6^z^ = 
0.002, the diffusion timescale is around three orbits. The diffusion timescale for radial 
diffusion is completely equivalent to equation (29). 

For a linear gravity field, DuUemond & Dominik (2004) derive a diffusion timescale 
similar to equation (29). On the other hand, Dubrulle et al. (1995) state a diffusion timescale 
of l/(i7oTf). This expression is actually a gravitational settling timescale that determines 
the amount of time it takes to increase the dust density in the mid-plane significantly due 
to gravity. Since the diffusion equilibrium sets in at very modest mid-plane overdensities for 
small dust grains, the timescale for such grains to reach diffusion equilibrium is much shorter 
than the gravitational settling timescale. 
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In Fig. 5 we plot the evolution of the logarithmic dust density averaged over the x- and 
y-directions for the run 64ajz. Starting at a time of zero orbits, curves are shown at two 
orbits time separation up to a time of ten orbits. The timescale to reach diffusion equilibrium 
is evidently around a few orbits (the saturated state is shown in Fig. 6). This is in good 
agreement with the analytical estimates given above. 



6.2. Measured Turbulent Diffusion Coefficients 

We now turn to measuring the turbulent diffusion coefficient from the equilibrium con- 
figuration that is illustrated in Fig. 4. According to equations (21) and (26), the equilibrium 
logarithmic dust density should be a cosine function if diffusion is the proper description of 
the turbulent transport. As an example of how the vertical diffusion coefficient is measured, 
we show in Fig. 6 the logarithmic dust density averaged over the radial and toroidal direc- 
tions for the run 64a_z at a time of t = 38 orbits. Also shown is the minimum cosine fit 
(dotted line). The fit is excellent, and this shows that here the turbulent transport of the 
dust grains is well-described as diffusion. In Fig. 7 we plot, for the same run, the full time 
evolution of the amplitude of the best-fit cosine function and the quality of the fit, Q. The 
fit quality is defined as 

^~ EJ(lnn).,(z.)]^ • ^ ^ 

Here the sum is taken over the entire vertical direction and {■ ■ ■)xy is used to denote the 
average taken over the x- and y-directions. In Fig. 7 we also plot the turbulent diffusion 
coefficient derived from the amplitude using equation (21). Both the amplitude and the 
turbulent diffusion coefficient stay approximately constant in time. The fit quality fluctuates 
by more than an order of magnitude, but is generally very good {Q is less than 0.02 at all 
times). A similar behavior is found for all runs with short friction times, both for vertical 
and for radial gravity. 

We have also run simulations without the short friction time approximation. Here the 
friction time must be at least a few times longer than the computational time-step of the gas 
in order to resolve the frictional acceleration in our explicit numerical scheme, but shorter 
than an orbital period for the fluid approach to be valid. We shall refer to such values 
of the friction time as intermediate friction times. Simulations with a freely evolving dust 
velocity serve both the purpose of showing in how far the short friction time approximation 
is valid, and also how the turbulent diffusion coefficient behaves when the friction time 
becomes larger and acceleration effects come into play. Remember that the short friction 
time approximation assumes that the dust grains can always reach an equilibrium velocity 
in one computational time-step. Hence effects such as vortex trapping in turbulent eddies 
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are not possible in the short friction time approximation. The time evolution of cosine 
amplitude, fit quality and turbulent diffusion coefficient for the intermediate friction time 
run 64c_x: (with radial gravity and a friction time of i7oTf = 0.02 corresponding to dust 
grains with radii of a few centimeters) is shown in Fig. 8. Here the amplitude changes a lot 
with time, and the fit quality Q rises above 0.1 on several occasions. Apparently diffusion is 
not at all times a good description of the turbulent transport in this run, even though the 
grains are still relatively well-coupled to the gas. Nevertheless, the time averaged diffusion 
coefficient in Fig. 7 and Fig. 8 is approximately the same as for the small grains. 

The measured turbulent diffusion coefficients in the vertical and radial directions are 
shown in Table 2. In Fig. 9 we plot the diffusion coefficients together with the a- value 
based on the Reynolds stress, at, and the a- value based on the Maxwell stress, We 
include l-cx fluctuation intervals on all measurements. The radial diffusion is seen to be 
much stronger (around 70%) than the vertical diffusion. This is also to be expected from 
Fig. 1, since the root-mean-square of the vertical velocity component is smaller than for the 
horizontal components, so the velocity fluctuations in the radial direction are stronger than 
in the vertical direction. 

From the length of the fluctuation bars in Fig. 9, it is clear that the fluctuations in 
the turbulent diffusion coefficient are very small for the short friction time limit. Combined 
with the fact that the quality of the cosine fit is excellent, this means that the turbulent 
transport in that case is well-described as diffusion. For intermediate friction times, with a 
freely evolving dust velocity, the fluctuation in the turbulent diffusion coefficient becomes 
larger, especially in the radial direction. The average values of the diffusion coefficients 
nevertheless stay approximately constant both for short and intermediate friction times. 
This gives confidence in that the short friction time approximation is indeed valid for very 
small dust grains. 



6.2.1. Diffusion Prandtl Number 

It is of great interest to compare the measured diffusion coefficients with the turbulent 
viscosity, since a popular parametrization of turbulent diffusion is to set the diffusion coef- 
ficient equal to the turbulent viscosity coefficient. It is seen from Table 2 that the vertical 
diffusion coefficient is generally around a factor of three to four times the non-magnetic a- 
value, but the value is comparable to the magnetic a-value. The radial diffusion coefficient 
is slightly higher than the total turbulent a-value. 

We quantify the difference between the measured turbulent diffusion coefficients and the 
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turbulent viscosity through the diffusion Prandtl number Pd. This is defined as the ratio 
between the turbulent viscosity and the turbulent diffusion coefficient as 

Pd=-^. (31) 

For anisotropic turbulence, the Prandtl number depends on the direction. Unfortunately 
there is no way to estimate the turbulent viscosity in the vertical direction, as there is no 
background shear and thus no flux of angular momentum vertically, so we shall use the value 
for the radial turbulent viscosity even for the vertical Prandtl number. 

The measured Prandtl numbers are shown in the last two columns of Table 2. The 
vertical diffusion Prandtl number is found to be above unity in the range Pd^ = 1.27 . . . 1.60, 
while the radial diffusion Prandtl number is below unity in the range Pd^^ = 0.79 . . . 0.90 and 
falling with increasing resolution. This is quite surprising as Prandtl numbers smaller than 
one could not be expected from standard diffusion theory. It is not possible to say whether 
the vertical Prandtl number would be similarly low if we had scaled with the proper vertical 
turbulent viscosity, because this quantity is, as mentioned above, not known. 



6.2.2. Dependence on Particle Size 

Much analytical work has been devoted to parametrizing the dependence of the diffusion 
coefficient on dust particle radius (Safronov 1969; Volk et al. 1980; Cuzzi et al. 1993; Dubrulle 
et al. 1995; Schrapler & Henning 2004; Reeks 2005). According to Schrapler & Henning 
(2004), ignoring the effect of the mean motion of the dust grains, the diffusion coefficient 
can be written as 

Here St is the Stokes number, and the factor 1/(1 + St) determines the variation of diffusion 
coefficient with particle radius. The Stokes number is defined as the ratio of the friction 
time to the turn-over time t^. of the largest eddies. Assuming that the rotation speed of the 
largest eddies as Ve = afcs and choosing q = 0.5, one can (following Schrapler & Henning 
2004) derive the expression 

Thus for the largest grains considered in this work, with l^o^f = 0.02, the expected change 
in diffusion coefficient due to particle size is around 1.5%. This is well below the fluctuation 
intervals in the measurements. Our results indeed confirm that there is no apparent size- 
dependence on the measured diffusion coefficient for our chosen grain size range. This 
also confirms the interpretation that the variation in the observed disc thickness at various 
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wavelengths is due to differential settling between particles of different sizes (e.g. DuUemond 
& Dominik 2004) and not due to a variation in the diffusion coefficient with particle size. 



6.3. Local Dust Density Enhancement 

To explore why the diffusion coefficient ffuctuates so much in the intermediate friction 
time runs, we plot in Fig. 10 the dust-to-gas ratio probability function for an intermediate 
friction time run with Qq^i = 0.02 (full line) and a short friction time run with QoTf = 2x10^'^ 
(dotted line) for a resolution of 64^ and no gravity. These two gravity-free runs are named 
64c_ng and 64a_ng, respectively. The probability of a grid point having a dust-to-gas ratio 
between and + Aea is 

, . A/(ed) 
Pie.) = , (34) 

where A/(ed) is the fraction of all grid points in the simulation box having a dust-to-gas 
ratio between and -|- Aej. We average over 10 orbits taken equidistantly between orbits 
10 and 100. According to Fig. 10 the probability of finding grid points with very high or very 
low dust-to-gas ratios is much higher in the intermediate friction time run than in the short 
friction time run. The dust-to-gas ratio in the short friction time run is extremely peaked 
around ed = cq, thus only the bottom part of the curve could be shown in Fig. 10. The full 
curve is shown in Fig. 11. 

There are two potential sources for the high dust-to-gas ratio contrast that is seen in 
the intermediate friction time: trapping of dust grains in turbulent vortices or trapping in 
regions of high pressure by pressure gradient trapping, as mentioned in the introduction. 
The latter effect can work also in the short friction time approximation, the first can not. 
According to equation (11), the terminal velocity of small dust grains climbing up the local 
pressure gradient is (we ignore gas velocity and set external gravity to zero) 

w = Ti[p-\VP - J X B)]. (35) 

The evolution of the dust number density of a fluid element is controlled by the continuity 
equation 

D In , , 

— = -V.«,. (36) 

where D/Dt = d/dt + [w ■ V) is the advective derivative of the flow. Combining equation 
(36) with equation (35) shows that dust should concentrate in regions where V ■ [p~'^(VP — 
J X B)] = V ■ -F < and be removed from regions where the divergence is negative. We 
examine whether this is the case in the two bottom panels of Fig. 12. Here the average 
dust-to-gas ratio (including 1-cr fluctuation intervals) is shown for bins in V ■ -F. The left 
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panel is for the short friction time run 64a_ng while the right panel is for the intermediate 
friction time run 64c_ng. For the intermediate friction time run, there is evidently some 
correlation between a positive divergence and a low dust-to-gas ratio and vice versa, but the 
correlation is not very strong. 

Vortex trapping is another potential source of the dust-to-gas ratio contrast (Barge & 
Sommeria 1995). It can be very powerful when i7oTf is close to unity. The delayed accel- 
eration of a dust grain entering a turbulent gas eddy causes the Coriolis force to dominate 
completely over the centrifugal force of the eddy. The effect of vortex trapping can be seen 
from the vorticity a; = V x w of the flow. Cyclonic vortices (with positive Uz) have an 
outwards directed Coriolis force relative to the center of motion and can expel dust grains. 
Anticyclonic vortices (with negative Uz) have a Coriolis force that points inwards. Such 
vortices can trap dust grains. As an illustration of the trapping of dust grains in turbulent 
features we show in Fig. 13 contour plots of Uz and ed in an arbitrarily chosen x-y-plane. The 
vorticity contours show patches of positive and negative vorticity. The correlation between 
negative vorticity and high dust-to-gas ratio (and vice versa) is clearly seen in many places. 
However, it is not a perfect 1:1 fit, as can also be expected in a dynamical system that 
is changing all the time. All concentrations are only surviving as long as a vortex exists. 
Turbulent eddies have a lifetime comparable to the shear time of the system, i.e. the orbital 
period. 

It is easier to see the correlation between vertical vorticity and dust-to-gas ratio in 
Fig. 12. Here the three top rows show the correlation between dust-to-gas ratio and the 
three directional components of the vorticity. There is a strong correlation with vertical 
vorticity component Uz for the intermediate friction time run. This is exactly as expected in 
case vortex trapping and expelling is the source of the number density contrast. A vertical 
vorticity can however also be caused by a non-rotating flow, e.g. if the gas-flow is hyper- 
Keplerian with a shear velocity that is linear with the radial coordinate Uy oc x. Such a 
profile can be caused by a radial bump in the gas density. Here dust-trapping would be due 
to pressure gradient trapping and not due to vortex trapping. 

A better test of vortex trapping than vertical vorticity can be devised by taking a closer 
look at the trapping mechanism (see e.g. Johansen et al. 2004). An anticyclonic vortex is in 
equilibrium because there is a resulting force on the gas particles pointing towards the center 
of rotation. This resulting force is a vector sum of the Coriolis force, the pressure gradient 
force and the Lorentz force, and it works as a centripetal force that supplies just the right 
amount of force necessary to orbit the center of rotation. In the fluid equations, the resulting 
centripetal force is balanced by the additional advection term that keeps the velocity field 
unchanged, even though the fiuid elements themselves experience an acceleration towards 
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the center of rotation. Thus for anticyclonic vortices, the advection vector —{u- V)w points 
away from the center of rotation, while the Coriohs force f{u) points towards the center 
of rotation, which is exactly in the opposite direction. The occurrence of the Coriolis force 
pointing in the opposite direction of the advection vector is a sufficient condition for having 
an anticyclonic vortex and thus vortex trapping. For a cyclonic vortex both the Coriolis 
force and the advection vector point away from the center of rotation. Defining the vortex 
parameter = [—{u- V)it] ■ f{u), we can now recognize cyclones by a positive value of iP^ 
and anticyclones by a negative value of \P. If dust grains are affected by vortex trapping, 
then there should be an anticorrelation between \P and the dust-to-gas ratio at the locations 
of cyclones and anticyclones. We examine this in Fig. 14. It is seen that the anticorrelation 
between ^ and dust-to-gas ratio is significant. This allows us to conclude that the large 
fluctuations in dust-to-gas ratio for the intermediate friction time runs is caused by trapping 
in turbulent eddies. 

Curiously there is also a significant correlation between any non-zero toroidal vorticity 
component Uy and a low dust-to-gas ratio for the intermediate friction time run in Fig. 12. 
This may be related to dust grains being expelled from eddies with a rotation axis parallel 
to the mid-plane (in the absence of gravity; when vertical gravity is included, particles can 
become suspended in such eddies, see Klahr & Henning 1997; Pasquero et al. 2003). However, 
there is no similar correlation with the radial component of vorticity Ux, probably because 
the shear wipes out any depletions/concentrations on a very short timescale. 

A similar search for concentrations of dust grains in MRI turbulence was performed by 
Hodgson & Brandenburg (1998). They find that for a frozen gas velocity field, intermediate 
friction time dust grains do indeed concentrate in the turbulent gas structures, but they 
attribute this effect to dust grains concentrating where the gas velocity field is converging 
rather than to vortex trapping. For an evolving gas velocity field, they find no concentration 
of dust. It is not clear why our results differ from these results. However, Hodgson & 
Brandenburg (1998) focus on concentrations of dust particles in the vertical plane, while 
in the current work dust concentrations are most pronounced in rotating structures in the 
horizontal plane. 

7. CONCLUSIONS 

The transport properties of dust grains in a turbulent accretion disc is of interest for 
many aspects of protoplanetary disc modeling and planet formation scenarios. In this paper, 
we have measured the turbulent diffusion coefficient of dust grains embedded in ideal MHD 
magnetorotational turbulence directly from numerical simulations. The choice of magne- 
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torotational turbulence was made because there is a growing realization that the magnetoro- 
tational instability can work at least in some parts of protoplanetary discs, even where the 
ionization fraction may be surprisingly low. It is also routinely produced in shearing box 
simulations, so it is a very accessible form of turbulence. Thus, by the use of MRI, we have 
a natural source of turbulence, whereas the current only other alternative for similar studies 
would be the use of driven turbulence in a box. The use of the ideal MHD equations can 
only be justified as a first approach to calculate the turbulent transport properties of dust 
grains. Further studies of non-ideal MHD should be made to clarify the transport properties 
of grains deeply embedded in the disc where the ionization fraction is low and where one is 
confronted with a "dead zone" around the mid-plane of the disc. 

As a numerical solver we have used the Pencil Code. This finite difference code solves 
the non-conservative form of the dynamical equations. It is special compared to other codes 
in that it uses sixth order derivatives in space. The numerical scheme of the Pencil Code 
was stabilized using hyperdiffusivity terms in all the dynamical equations. The effect of 
hyperdiffusivity is to affect the large scale motion as little as possible, while at the same 
time quenching unstable modes at the smallest scales of the box. By varying the size of 
the artificial diffusion coefficient, we have found the direct influence of artificial diffusion 
on the measured turbulent diffusion coefficient to be negligible, most likely due to the fact 
that mass diffusion is primarily contributed by the fast and far moving large scales of the 
turbulence, and these are as mentioned affected only very little by hyperdiffusivity. From 
this perspective, hyperdiffusivity seems to be a tool that is well suited for measurements of 
turbulent transport properties. 

Since we have only considered dust grains of sizes much less than one meter, the dust 
grains could be treated as a fluid interacting with the gas through a drag force. For the 
tiniest dust grains, where the friction time is much shorter than the computational time- 
step, we have used an algebraic equation to obtain the dust velocity at each time-step. This 
short friction time approximation incorporates the tendency of dust grains to move up the 
local pressure gradient of the gas, an effect which we have referred to as pressure gradient 
trapping. It can explain such phenomena as the settling of dust grains towards the mid-plane 
of a stratified disc and the radial drift of dust grains in discs with a radial pressure gradient. 
In the current work, we have also included the effects of magnetic pressure and tension in 
the short friction time approximation. For intermediate friction time dust grains, where the 
friction time is within a few orders of magnitude of the orbital period, we have integrated 
the dust equation of motion together with the other dynamical equations. Here acceleration 
effects are allowed, in the sense that dust grains are no longer assumed to instantaneously 
reach a terminal velocity where the drag force is balanced by the other forces affecting the 
dust. The fact that the time average of the measured diffusion coefficient was approximately 



-27- 



the same for tiny dust grains, using the short friction time approximation, and intermediate 
size dust grains, with a free evolution of the dust velocity, gives some credit to the validity 
of the short friction time approximation. 

We have chosen to measure the turbulent diffusion coefficient by forcing the dust grains 
to settle towards a mid-plane by an external force field. This settling was eventually balanced 
by turbulent diffusion away from the mid-plane. To deduce the value of the turbulent 
diffusion coefficient, the equilibrium dust density could then be compared with an analytical 
solution for a parametrized diffusion coefficient. The method works not only for the vertical 
direction, but also for the radial direction, so that we have been able to measure both the 
vertical and the radial turbulent diffusion coefficients. 

For the short friction time runs, the equilibrium dust number density was excellently 
fitted with the expected analytical solution. That means that the turbulent transport of small 
dust grains is well-described as diffusion. For intermediate friction times, the equilibrium 
dust number density was much more erratic, especially in the radial direction, and did not 
always give a good fit. We also found that the dust-to-gas ratio probability distribution was 
much wider than in the short friction time runs. To explore the reason for the large spread 
in dust-to-gas ratio, we have examined correlations between different parameters of the gas 
and the dust-to-gas ratio. A strong correlation between vertical vorticity component and 
dust-to-gas ratio was found. Based on this and an equally strong correlation between the 
sign of the vortex parameter and the dust-to-gas ratio, we conclude that the spread in dust- 
to-gas ratio, and thus the fluctuations in the diffusion coefficient, is due to vortex trapping in 
turbulent eddies (Barge & Sommeria 1995). Some weaker indications that pressure gradient 
trapping is taking place were also found, but similar to the results of Johansen et al. (2004), 
the over all dominant trapping mechanism is found to be vortex trapping. The dust-trapping 
that is seen in the current work happens for relatively well-coupled particles with a friction 
time on the order of a few percent times the shear time f^Q^. One can speculate that for 
larger particles, dust-trapping mechanisms will be so efficient that the diffusion picture of 
turbulent transport will no longer be valid, but further investigations into the transport of 
larger particles will have to examine this. Such an investigation would have to incorporate 
dust grains as particles moving on top of the gas fluid, since the fluid description of dust 
grains is no longer valid when the mean free path becomes larger than the scale height of 
the disc. 

In the vertical direction the turbulent diffusion coefficient was measured to be smaller 
than the total turbulent viscosity and have a diffusion Prandtl number of approximately 
Pd^ = 1.5. The diffusion coefficient is still considerably larger than the non- magnetic tur- 
bulent viscosity alone. The measured radial turbulent diffusion coefficient turned out to be 
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almost twice as large as the vertical diffusion coefficient. It is systematically larger than 
the total turbulent viscosity, i.e. the sum of the non-magnetic and the magnetic turbulent 
viscosity, with a diffusion Prandtl number of around Pd^ = 0.85. The value of the radial 
diffusion Prandtl number was found to be falling with increasing resolution. Future simula- 
tions should try to find convergence for the diffusion Prandtl numbers, but this is beyond 
the scope of the present work. 

The anisotropy between the vertical and the radial directions should be taken into 
account for studies of planetesimal formation which invoke a gravitational instability in the 
dust sublayer. Here the onset of a gravitational instability in the vertically settled dust 
layer depends strongly on the effect of vertical diffusion. The amount of anisotropy can 
be expected to increase if the effect of vertical gravity and stratification is included (S. 
Fromang, personal communication), since then the buoyancy of the gas would decrease the 
vertical velocity fluctuations. 

We want to stress that even though we find a radial diffusion Prandtl number of less than 
unity, the disc can still be assumed to radially transport dust grains and angular momentum 
about equally well. This is important for the modeling of radial mixing of dust grains and 
chemical species, a task which is becoming ever more relevant as more observations of the 
radial distribution of dust grains and molecules in protoplanetary discs become available. 
The result is in agreement with Yousef et al. (2003), who find for simulations of forced MHD 
turbulence a turbulent magnetic Prandtl number of unity. For dust grains, the equality 
between the radial turbulent diffusion coefficient and the turbulent viscosity is surprising, 
when considering that most of the angular momentum is transported by magnetic Maxwell 
stresses, while the dust grains have no coupling with magnetic fields at all. Following the 
argument of Tennekes & Lumley (1972), both angular momentum transport by Reynolds 
stresses and radial diffusion depend on the radial velocity fluctuations, so one would expect 
the non-magnetic a-value and the diffusion coefficient to be similar. The actual cause of the 
measured mismatch between turbulent diffusion and non-magnetic turbulent viscosity would 
seem to need more discussion in the future. 
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A. HYPERDIFFUSIVITY 



In this appendix we discuss the use of hyperdiffusivity and present the hyperversions of 
viscosity, mass diffusion and resistivity that we are applying in this work. 

Because the Pencil Code is a finite-difference code, artificial diffusivity terms are needed 
in all dynamical equations to stabilize the numerical scheme. For this purpose, we use sixth 
order hyperdiffusivity terms which affect mainly high wave numbers, the smallest scales in 
the simulation, and preserve the energy at low wave numbers. Hyperviscosity and magnetic 
hyperdiffusivity have been used extensively to study the properties of forced magnetohy- 
drodynamical turbulence (e.g. Brandenburg & Sarson 2002, and references therein). The 
prospect is to affect large scales as little as possible by dissipation, thus widening the inertial 
range beyond what can be achieved with a regular viscosity operator. 

Possible side effects of using hyperviscosity and magnetic hyperdiffusivity is to increase 
the bottleneck effect (a physical effect in turbulence where energy piles up around the dissi- 
pative scale, see e.g. Miiller & Biskamp 2000) and to cause the dynamo-generated magnetic 
field in helical flows to saturate at a higher level than what is seen when using a regular 
viscosity operator (Brandenburg & Sarson 2002). Nevertheless, for forced non- magnetic tur- 
bulence Haugen & Brandenburg (2004) show that the shape of the inertial range for runs with 
hyperviscosity is very similar to the shape for higher resolution runs without hyperviscosity. 

For the current work we define a momentum-conserving hyperviscosity function f^^ as 



Here S^^ is a simplified Ith order rate-of-strain tensor defined as 



For computational simplicity we consider the dynamical viscosity = to be constant. 
Then the hyperviscosity function takes the appearance 



where V^"^ = V^™' + V^™' -|- V^'" is a high order differential operator that reduces to a 
Laplacian for m = 1. For the purpose of stabilizing the numerical scheme we adopt a sixth 
order hyperviscosity by setting m = 3 in equation (A3). The hyperviscosity function 
then appears as 




(Al) 




(A3) 




(A4) 
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For the artificial mass diffusion term we define the hyperdiffusion function /d as 

/D(p) = (-ir-^/^™V2>, (A5) 

where is a constant diffusion coefficient. Using /d as a diffusion term in the continuity 
equation conserves mass density. Again we adopt a hyperdiffusivity version with m = 3, 
leading to the expression 

/d(p)=/^3VV. (A6) 

The hyperresistivity function is defined as 

/,(A) = (-1)™"S„V2-A, (A7) 

where rjm is the magnetic diffusivity. As for viscosity and diffusion we use a hyperresistivity 
scheme with m = 3 in equation (A7). Then the resistivity function f comes out as 

/,(A)=r/3V6A. (A8) 

Using this function as a resistivity term in the induction equation conserves all components 
of the magnetic field B. 



B. RADIAL DIFFUSION EQUILIBRIUM 

In this appendix we derive the equilibrium dust density when the dust is exposed to a 
radial gravity. We define a sinusoidal gravity field similar to what was done in the ^-direction 

as 

gx = -9QSin{k^x) , (Bl) 

where kx = 27t/Lx is the radial wave number of the field. In the horizontal plane the Coriolis 
force connects the radial and toroidal motions, so that any velocity in one direction results 
in an acceleration in the other direction. If not for the damping effect of friction, dust grains 
starting with any non-zero velocity would be forced to move in epicyclic motion. Fortunately 
the drag force from the gas permits an equilibrium solution to the dust equation of motion 
equation (5). We solve for u = Wz = 0, Wx = Wx{x) and Wy = Wy{x) and get the two 
component equations of the equation of motion 

Ow 1 
= -Wx^rr^ + 2f2oWy Wx - go sm{kxx) , (B2) 

ox Tf 

Ow 1 1 

= -Wx^r^ - -f^oWx Wy. (B3) 

dx 2 T{ 
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Again the advection is ignored. This leads to an algebraic linear system of equations in Wx 

and w„ 








2f2oWy Wx - gosm{kxx) 

Tf 

1 ^ 1 



-f2oWx 



-w. 



y ' 



(B4) 
(B5) 



that has the solution 



—MO 



sm{kxx) 



^QTf go 



2{1+I22r2) 



At sin(A;^.x) 



Tigo sm{kxx) 
QoT^qq sin(/c^x) 



(B6) 



Here the approximate expressions are valid to first order in QoTf. The ratio between toroidal 
and radial velocity is \wy/wx\ = |^oTf, so the toroidally imposed velocity becomes unim- 
portant with sufficiently short friction time. For this form of velocity field, the equilibrium 
continuity equation takes the form 



- - — 

dx 



d_ 

dy 



dy 



(B7) 

By considering solutions to equation (B7) of the form n{x,y) = n{x), the ?/-derivative term 
of equation (B7) vanishes entirely. Then the equilibrium solution to the continuity equation, 
when assuming that the total fiux of number density radially through the box is zero, is 
simply 



In n = In ni + 



_Tj9o_ 

krrDx ^ 



cos{kxx) 



(B8) 



formally identical to the vertical case. 
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Fig. 1. — Evolution of various turbulence parameters for a 64^ run (left panels) and a 
128^ run (right panels). The top panels show the evolution of total kinetic energy and its 
directional components. The radial and toroidal directions have comparable values of kinetic 
energy, whereas the vertical direction has around a factor of two less. The magnetic energy 
(middle panels) is completely dominated by the toroidal magnetic field. The uy component 
of the Reynolds and Maxwell stresses (lower panels) is effectively a measure of the turbulent 
viscosity. The magnetic stresses are around four times higher than the kinetic stresses. 



- 36 - 




Fig. 2. — Fourier spectrum of the velocity components of the gas for a 64^ resolution run, 
averaged from 10 to 100 orbits. A Kolmogorov k~^/^ line is shown for reference. Both the 
radial and the toroidal components show a Kolmogorov-like behavior on large scales, whereas 
the vertical component is flatter. At small scales dissipation becomes important. The radial 
and vertical directions show a rise in power on the very smallest scales. 




Fig. 3. — Same as Fig. 2, but for a 128^ resolution run. The vertical velocity component 
still has a flatter slope in the inertial range than the two other components, but the surplus 
power at the very smallest scale is greatly diminished compared to the 64^ run. 
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Fig. 4. — Dust density contours at the sides of the simulation box for the short friction time 
run 128a_z. The radial direction is towards the right while the shearing direction is towards 
left. The dust is concentrated around the mid-plane due to a vertical gravity acting only on 
the dust. Turbulent transport alone prevents the further vertical settling of the dust layer. 
This configuration is statistically unchanged for at least one hundred orbits. 
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Fig. 5. — The logarithmic dust density of the run 64a_z, averaged over x and at different 
times. The curves are each separated by two orbits going from t = (full line) to t = 10 orbits 
(long-dashed line). The approach to equilibrium happens on a timescale of a few orbits, in 
good agreement with the analytical estimate of the diffusion timescale that is presented in 
the text. 
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Fig. 6. — The logarithmic dust density averaged over x and y as function of vertical height 
z (full line) and a cosine fit (dotted line). The cosine fit is in excellent agreement with the 
data. This shows that the turbulent transport is indeed well described as diffusion. Shown 
here is for the short friction time run 64a_z at a time of 26 orbits. The fit quality (defined 
in the text) is Q ^ 0.005. 
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Fig. 7. — Time evolution of the fitted cosine amplitude (dotted line), the fit quality (dash- 
dotted line) and the derived radial turbulent diffusion coefficient (full line) for the short 
friction time run 64a_x: with radial gravity. Both the fit amplitude, and hence the turbulent 
diffusion coefficient, are approximately constant in time, although small variations are seen. 
The fit quality is generally excellent, but it fluctuates with around an order of magnitude 
during the 100 orbits shown here. Compare with Fig. 8 which shows the evolution of the 
same variables for an intermediate friction time run. 
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Fig. 8. — Same as Fig. 7, but for the intermediate friction time run 64c_x. Obviously the 
cosine amphtude and the derived turbulent diffusion coefficient change much more violently 
with time. The fit quality is also a lot worse. The average diffusion coefficient is actually 
the same as for the short friction time run shown in Fig. 7, but the poor fit quality here 
means that the diffusion description of turbulent transport is not as good as it is in the short 
friction time runs. 
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Fig. 9. — The measured turbulent diffusion coefficient as a function of i7oTf for a resolution 
of 64^ (left panels) and 128^ (right panels). The vertical diffusion coefficient is shown in the 
top panels, while the radial diffusion coefficient is shown in the bottom panels. For reference 
the turbulent a-values based on both the Reynolds stress and the Maxwell stress are shown 
including their 1-a fluctuation intervals. The radial diffusion coefficient is comparable to the 
sum of the turbulent a- values and is around 70% higher than the vertical diffusion coefficient. 
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Fig. 10. — Dust-to-gas ratio probability distribution function for runs with O^ri = 0.02 (full 
line) and OoTf = 2 x 10"'' (dotted line) without gravity on the dust. For the intermediate 
friction run, there is a much higher probability for very low or very high dust-to-gas ratios, 
compared to the short friction time run where the dust-to-gas ratio is sharply peaked around 
ed = Cq. The full probability curve for the short friction time run is shown in Fig. 11. 
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Fig. 11. — The full dust-to-gas ratio probability distribution function for a short friction 
time run with f^oTf = 2 x 10~^ and no gravity. Because the peak is so sharp compared to 
the intermediate friction time run, only the lower part is shown in Fig. 10. 
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Fig. 12. — The dust-to-gas ratio in bins of vorticity components (first three rows) and 
divergence of pressure gradient flux (last row). The large dot shows the average dust-to-gas 
ratio in the bin, while the bars represent the fluctuation interval. The clearest correlation is 
between vertical vorticity and dust-to-gas ratio for the intermediate friction time run. This 
may be due to vortex trapping as explained in the text. 
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Fig. 13. — Contour plots of vertical component of vorticity (left panel) and dust-to-gas ratio 
(right panel) in an arbitrary 2;-plane for the intermediate friction time run 64c_ng. There is 
a tendency for positive vorticity (light regions) to correspond to low dust-to-gas ratio (dark 
regions) and vice versa. This indicates that dust grains are being trapped in turbulent eddies 
by the vortex trapping mechanism. The dotted lines are reference lines to make comparison 
between the two plots easier. 



- 46 - 




1.02 



1.01 - 



1.00 



0.99 - 



0.98 

-0.0010 -0.0005 0.0000 0.0005 0.0010 -0.0005 0.0000 



0.0005 0.0010 



Fig. 14. — Plot of dust-to-gas ratio in bins of vortex parameter ^ = [—{u ■ 'V)u] ■ f{u). 
Anticyclonic vortices have a negative value of \l/ , whereas for cyclonic vortices ip' is positive. 
For the intermediate friction time run, there is a clear anticorrelation between vortex param- 
eter and dust-to-gas ratio. This is an indication that dust is being trapped in anticyclonic 
vortices. 
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Table 1. Run parameters 



Run 


Res 


(0) 

9z 


yx 








a,/m 




SFTA 


Ai3 




■ V3 = 






643 


1000 





2 


X 10' 


-7 


2 


X 


10' 


'7 


Yes 


1 


.3 


X 


10- 


11 


64b^ 


643 


10 





2 


X 10' 


-5 


2 


X 


10' 


-5 


Yes 


1 


.3 


X 


10' 


11 


64c_z 


643 


0.01 







0.02 






0. 


,02 




No 


1, 


.3 


X 


10' 


11 


128aj 


1283 


1000 





2 


X 10' 


-7 


2 


X 


10' 


-7 


Yes 


1 


.3 


X 


10' 


-12 


128c_z 


1283 


0.01 







0.01 






0. 


,01 




No 


1, 


.3 


X 


10' 


-12 


64a_x; 


643 





1000 


2 


X 10' 


-7 


2 


X 


10' 


-7 


Yes 


1 


.3 


X 


10' 


-11 


64b^ 


643 





10 


2 


X 10' 


-5 


2 


X 


10' 


-5 


Yes 


1 


.3 


X 


10' 


-11 


64c_x 


643 





0.01 




0.02 






0, 


,02 




No 


1, 


.3 


X 


10' 


-11 


128a^ 


128^ 





1000 


2 


X 10' 


-7 


2 


X 


10' 


-7 


Yes 


1 


.3 


X 


10' 


-12 


128c^ 


1283 





0.01 




0.01 






0, 


,01 




No 


1, 


.3 


X 


10' 


-12 


64a_ng 


643 








2 


X 10' 


-7 


2 


X 


10' 


-7 


Yes 


1 


.3 


X 


10' 


-11 


64c_ng 


643 










0.02 






0. 


,02 




No 


1, 


.3 


X 


10' 


-11 



Note. — The first column gives the name of the run, the second the reso- 
lution, the third and fourth the vertical and radial gravity strength, the fifth 
column the friction time, the sixth column the corresponding grain radius in 
a typical solar nebula at ro = 5 AU, the seventh column whether we used the 
short friction time approximation or not, and the eighth column the value of 
the artificial viscosity /xs, magnetic diffusivity r/a and mass diffusion D^. 
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Table 2. Turbulent viscosities and turbulent diffusion coefficients 



Run 


at/10-=^ 








Pd. 


Pd. 


64a_z 


0.34 ±0.07 


1.52 ±0.27 


1.18±0.11 




1.58 




64b_z 


0.34 ±0.07 


1.52 ±0.27 


1.18 ±0.11 




1.58 




64c_z 


0.33 ±0.08 


1.47 ±0.29 


1.12 ±0.14 




1.60 




64a_x; 


0.34 ±0.07 


1.52 ±0.27 




2.07 ±0.28 




0.90 


64b^ 


0.34 ±0.07 


1.52 ±0.27 




2.07 ±0.28 




0.90 


64c_x 


0.33 ±0.08 


1.47 ±0.29 




2.12 ±0.75 




0.85 


128a_z 


0.19 ±0.03 


0.85 ±0.12 


0.82 ±0.10 




1.27 




128c_z 


0.19 ±0.04 


0.85 ±0.19 


0.79 ±0.13 




1.31 




128a^ 


0.16 ±0.02 


0.75 ±0.10 




1.15 ±0.14 




0.79 


128CJC 


0.18 ±0.03 


0.83 ±0.12 




1.27 ±0.30 




0.79 



Note. — The first column gives the name of the run. The second and 
third columns show the turbulent a-values based on Reynolds and Maxwell 
stresses, respectively. Since there is no back-reaction from the dust on the gas, 
these values are only affected by the dust through the dust's contribution to 
the time-step. The next two columns show the measured turbulent diffusion 
coefficients, and in the last two columns we write the vertical and radial 
turbulent diffusion Prandtl numbers. 



